function [out] = Gaussian(in, fc)
% Gaussian filter from 
% Colquhoun and Sigworth (Pratical Analysis of Records, Chapter 19 in
% Sakmann and Neher, Single Channel Recording).
% They claim that this is the optimal filter for the first
% pass of analysis.
% in is the input array
% fc is the -3db cutoff frequency, in units of sampling frequency
% 12/29/99 P. Manis
%
sigma = 0.132505/fc;
if(sigma<0.62)
	a(2)=sigma*sigma/2.0; % narrow impulse response - only 3 terms used
	a(1)=1.0-2*a(2);
	nc=1;
else 
	% sigma > 0.62
	% standard gaussian coefficients.
	% nc is the number of coeffs not counting the central one a(0)
	%
	nc = floor(4.0*sigma);
	if(nc > 53)
		nc = 53;
	end;
	b = -0.5/(sigma*sigma);
	a(1)=1.0;
	sum=0.5;
	
	for i=1:nc
		t=exp(i*i*b);
		a(i+1)=t;
		sum=sum+t;
	end;
	% normalize the coefficients
	%sum=sum*2;
	for i=1:nc+1
		a(i)=a(i)/sum;
	end;
end;
%
%
[nr, np] = size(in);
for l=1:nr % across all records
	out(l,:)=conv(in(l,:),a); % do filter by convolution

% actual filter operates here. This implementation is EXACTLY like sigworth's
% (except for handling the two dimensional array directly)
%	for i=1:np
%		jl=i-nc;
%		if(jl < 1)
%			jl = 1;
%		end;
%		ju=i+nc;
%		if(ju > np)
%			ju = np;
%		end;
%		sum = 0;
%		for j=jl:ju
%			k = abs(j-i)+1;
%			sum = sum+in(l,j)*a(k);
%		end;
%		out(l,i)=sum(in(l,[jl:ju]).*a[jl-i+1:ju-i+1])
%		out(l,i)=sum;
%	end;
end;

		
